##################################################################
# /* Author: Gautam Nair and Nicholas Sambanis */
# /* Violence Exposure and Ethnic Identification: Evidence from Kashmir */
# /* Appendix Table 7-10: Correlation Between Self Reported Violence Exposure and Identification  */

##################################################################

# change working directory here
# setwd("")

##################################################################
# loading packages
##################################################################
library(Hmisc)
library(foreign)
library(sandwich)
library(lmtest)
library(numDeriv)
library(stargazer)
library(ggplot2)
library(plyr)
library(gridExtra)
library(dplyr)
library(scales)

##################################################################

rm(list=ls())

library(stargazer)

data.working <- read.csv("nei_kashmir_replication_dataset.csv")

dep.var <- c("identity_rank_indian_rev", 
"identity_choice_indian_both", 
"kashmir_status_india",  
"ind_pak_prefer_ind",
"bonus_donation_india", 
"ic_index_id_no_ind_pak",
"ic_index_id_ind_pak",
"protest_peaceful_endorse",
"protest_violent_endorse"
)

ind.violence <- c("violence_exp_index_dummy",
				"group2_violence",
				"group3_growth",
				"group4_int_inst",
                  "group5_geography")

ind.correlates <- c("age_18_27",               
                    "age_28_49",               
                    "gender_female",
                    "education_class12",
                    "religion_sunni",
                    "religion_shia",
                    "main_earner",
                    "index_income_wealth1",
                    "occupation_govt")

ind.districts <- c("constituency_1",  
                   "constituency_2", 
                   "constituency_3", 
                   "constituency_4", 
                   "constituency_5", 
                   "constituency_6", 
                   "constituency_7", 
                   "constituency_8", 
                   "constituency_9", 
                   "constituency_10", 
                   "constituency_11", 
                   "constituency_12", 
                   "constituency_13", 
                   "constituency_14", 
                   "constituency_15", 
                   "constituency_16", 
                   "constituency_17", 
                   "constituency_18", 
                   "constituency_19", 
                   "constituency_20", 
                   "constituency_21", 
                   "constituency_22")
# constituency X is excluded

ind.1 <- ind.violence
ind.2 <- c(ind.violence, ind.correlates)
ind.3 <- c(ind.violence, ind.correlates, ind.districts)

reg.models <-vector("list", length(dep.var)) 
se.models <-vector("list", length(dep.var))
matrix <- vector("list", length(dep.var))
x=0

for(i in 1:length(dep.var)){
  
  my.formula <-paste(dep.var[i],'~', paste(ind.1, collapse= ' + '))
  temp.model <- lm(my.formula, data=data.working)
  temp.se <- coeftest(temp.model, vcov=vcovHC(temp.model,type="HC2"))[,2]
  x = x+1
  se.models[[x]] <- temp.se
  reg.models[[x]] <- temp.model
  
  my.formula <-paste(dep.var[i],'~', paste(ind.2, collapse= ' + '))
  temp.model <- lm(my.formula, data=data.working)
  temp.se <- coeftest(temp.model, vcov=vcovHC(temp.model,type="HC2"))[,2]
  x = x+1
  se.models[[x]] <- temp.se
  reg.models[[x]] <- temp.model	
  
  my.formula <-paste(dep.var[i],'~', paste(ind.3, collapse= ' + '))
  temp.model <- lm(my.formula, data=data.working)
  temp.se <- coeftest(temp.model, vcov=vcovHC(temp.model,type="HC2"))[,2]
  x = x+1
  se.models[[x]] <- temp.se
  reg.models[[x]] <- temp.model	
  
}


varlabels <- c("Violence Exposure",
               "Age: 18-27", 
               "Age: 28-49", 
               "Gender: Female",
               "$>=$High School",
               "Religion: Sunni",
               "Religion: Shia",
               "Main Earner",
               "Wealth Index",
               "Occupation Govt")

dropvars <- c("group2_violence","group3_growth","group4_int_inst","group5_geography", ind.districts)

spectitle <- c("Observed Correlation between Self-Reported Violence Exposure and Identification: Part 1/4")
outputfile <- c("tf_t_07_violence_identification_correlation_1")
cutoffs.star <- c(0.05, 0.01)
char.star = c("**", "***")
clabels <- c("Ranked Indian Identity (1-4)", "Identity: Indian/Indian+Kashmiri (1/0)")
cseparate <- c(3,3)
temptype <- c("latex", "text")
tempext <- c(".tex", ".txt")

for(q in 1:length(temptype)){
	temp.output <- paste(outputfile, tempext[q], sep="")
stargazer(reg.models[[1]], reg.models[[2]], reg.models[[3]], 
	reg.models[[4]], reg.models[[5]], reg.models[[6]], 
          se=list(se.models[[1]], se.models[[2]], se.models[[3]], 
                  se.models[[4]], se.models[[5]], se.models[[6]]),                
          title= spectitle,
          out=temp.output,
          no.space=TRUE, 
          model.numbers= TRUE,
          notes.append=FALSE,
          align=FALSE,
          covariate.labels=varlabels,
          header= FALSE,
          font.size="scriptsize",
          model.names=FALSE,
          digits=2,
          omit.stat = c("rsq", "f","adj.rsq"),
          dep.var.caption  = "",
          dep.var.labels.include=FALSE,
          df=FALSE,
          omit.table.layout="n",
          star.cutoffs = cutoffs.star,
          star.char = char.star,
          order=ind.3,
          omit=dropvars,
          report="vc*s",
          column.sep.width="0.1pt",
          column.labels = clabels,
          column.separate = cseparate,
          type=temptype[q],
           add.lines = list(c("Individual Covariates", "No", "Yes", "Yes", "No", "Yes", "Yes"),
                           c("Constituency Fixed Effects", "No", "No", "Yes", "No", "No", "Yes"))
)
}

spectitle <- c("Observed Correlation between Self-Reported Violence Exposure and Identification: Part 2/4")
outputfile <- c("tf_t_08_violence_identification_correlation_2")
cutoffs.star <- c(0.05, 0.01)
char.star = c("**", "***")
clabels <- c("Kashmir Status: Remain in India (1/0)", "India/Pakistan: Prefer India (1/0)",  "Donation to Indian NGO (0-100)")
cseparate <- c(3,3,3)
temptype <- c("latex", "text")
tempext <- c(".tex", ".txt")

for(q in 1:length(temptype)){
	temp.output <- paste(outputfile, tempext[q], sep="")
stargazer(reg.models[[7]], reg.models[[8]], reg.models[[9]],
reg.models[[10]], reg.models[[11]], reg.models[[12]], 
	reg.models[[13]], reg.models[[14]], reg.models[[15]],  
          se=list(se.models[[7]], se.models[[8]], se.models[[9]],
          se.models[[10]], se.models[[11]], se.models[[12]], 
                  se.models[[13]], se.models[[14]], se.models[[15]]),                
          title= spectitle,
          out=temp.output,
          no.space=TRUE, 
          model.numbers= TRUE,
          notes.append=FALSE,
          align=FALSE,
          covariate.labels=varlabels,
          header= FALSE,
          font.size="scriptsize",
          model.names=FALSE,
          digits=2,
          omit.stat = c("rsq", "f","adj.rsq"),
          dep.var.caption  = "",
          dep.var.labels.include=FALSE,
          df=FALSE,
          omit.table.layout="n",
          star.cutoffs = cutoffs.star,
          star.char = char.star,
          order=ind.3,
          omit=dropvars,
          report="vc*s",
          column.sep.width="0.1pt",
          column.labels = clabels,
          column.separate = cseparate,
          type=temptype[q],
		add.lines = list(c("Individual Covariates", "No", "Yes", "Yes", "No", "Yes", "Yes", "No", "Yes", "Yes"),
		 c("Constituency Fixed Effects", "No", "No", "Yes", "No", "No", "Yes", "No", "No", "Yes"))
)
}

spectitle <- c("Observed Correlation between Self-Reported Violence Exposure and Identification: Part 3/4")
outputfile <- c("tf_t_09_violence_identification_correlation_3")
cutoffs.star <- c(0.05, 0.01)
char.star = c("**", "***")
clabels <- c("National Identification Index I", "National Identification Index II")
cseparate <- c(3,3)
temptype <- c("latex", "text")
tempext <- c(".tex", ".txt")

for(q in 1:length(temptype)){
	temp.output <- paste(outputfile, tempext[q], sep="")
stargazer(reg.models[[16]], reg.models[[17]], reg.models[[18]], 
	reg.models[[19]], reg.models[[20]], reg.models[[21]], 
          se=list(se.models[[16]], se.models[[17]], se.models[[18]], 
                  se.models[[19]], se.models[[20]], se.models[[21]]),                
          title= spectitle,
          out=temp.output,
          no.space=TRUE, 
          model.numbers= TRUE,
          notes.append=FALSE,
          align=FALSE,
          covariate.labels=varlabels,
          header= FALSE,
          font.size="scriptsize",
          model.names=FALSE,
          digits=2,
          omit.stat = c("rsq", "f","adj.rsq"),
          dep.var.caption  = "",
          dep.var.labels.include=FALSE,
          df=FALSE,
          omit.table.layout="n",
          star.cutoffs = cutoffs.star,
          star.char = char.star,
          order=ind.3,
          omit=dropvars,
          report="vc*s",
          column.sep.width="0.1pt",
          column.labels = clabels,
          column.separate = cseparate,
          type=temptype[q],
		 add.lines = list(c("Individual Covariates", "No", "Yes", "Yes", "No", "Yes", "Yes"),
		c("Constituency Fixed Effects", "No", "No", "Yes", "No", "No", "Yes"))
)
}

spectitle <- c("Observed Correlation between Self-Reported Violence Exposure and Identification: Part 4/4")
outputfile <- c("tf_t_10_violence_identification_correlation_4")
cutoffs.star <- c(0.05, 0.01)
char.star = c("**", "***")
clabels <- c("Peaceful Protest Participation (1-5)",  "Violent Protest Endorsement (1-5)")
cseparate <- c(3,3)
temptype <- c("latex", "text")
tempext <- c(".tex", ".txt")

for(q in 1:length(temptype)){
	temp.output <- paste(outputfile, tempext[q], sep="")
stargazer(reg.models[[22]], reg.models[[23]], reg.models[[24]], 
	reg.models[[25]], reg.models[[26]], reg.models[[27]], 
          se=list(se.models[[22]], se.models[[23]], se.models[[24]], 
                  se.models[[25]], se.models[[26]], se.models[[27]]),                
          title= spectitle,
          out=temp.output,
          no.space=TRUE, 
          model.numbers= TRUE,
          notes.append=FALSE,
          align=FALSE,
          covariate.labels=varlabels,
          header= FALSE,
          font.size="scriptsize",
          model.names=FALSE,
          digits=2,
          omit.stat = c("rsq", "f","adj.rsq"),
          dep.var.caption  = "",
          dep.var.labels.include=FALSE,
          df=FALSE,
          omit.table.layout="n",
          star.cutoffs = cutoffs.star,
          star.char = char.star,
          order=ind.3,
          omit=dropvars,
          report="vc*s",
          column.sep.width="0.1pt",
          column.labels = clabels,
          column.separate = cseparate,
          type=temptype[q],
		add.lines = list(c("Individual Covariates", "No", "Yes", "Yes", "No", "Yes", "Yes"),
	 c("Constituency Fixed Effects", "No", "No", "Yes", "No", "No", "Yes"))
)
}





